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§13.2 GENERALIZED INTERPOLATION 


§13.1. Introduction 

This Chapter develops special one-dimensional elements, such as thick beams, arches and beams on elastic 
foundations, that require mathematical and modeling resources beyond those presented in Chapters 11-12. 
The techniques used are less elementary, 1 and may be found in books on Advanced Mechanics of Materials, 
e.g. [97,108]. Readers are expected to be familiar with ordinary differential equations and energy methods. 
The Chapter concludes with beam accuracy analysis based on the modified equation method. 

All of this Chapter material would be normally bypassed in an introductory finite element course. It is primarily 
provided for offerings at an intermediate level, for example a first graduate FEM course in Civil Engineering. 
Such courses may skip most of Part I as being undergraduate material. 

§13.2. Generalized Interpolation 

For derivation of special and C° beam elements it is convenient to use a transverse-displacement cubic inter¬ 
polation in which the nodal freedoms iq, t> 2 , 0\ and 0 2 are replaced by generalized coordinates c\ to c 4 : 

v($) = N c i c 1 + N C 2 C2 + N c 3 c 3 + N c 4 C4 = N c c. (13.1) 

Here N n (£) are generalized shape functions that satisfy the completeness requirement discussed in Chapter 19. 
N, is a 1 x4 matrix whereas c is a column 4-vector. Formula (13.1) is a generalized interpolation. It includes 
the Hermite interpolation (12.10-12.12) as an instance when c 1 = tq, c 2 = v[, c 3 = v 2 and c 4 = v' 2 . 

§13.2.1. Legendre Polynomials 

An obvious generalized interpolation is the ordinary cubic polynomial v(f) — <q + c 2 § + c 3 § 2 + c 4 § 3 , but this 
turns out not to be particularly useful. A more seminal expression is 

u(£) = L 1 <q + L 2 c 2 + L 3 c 3 + L 4 c 4 = L c, (13.2) 


where the L, are the first four Legendre polynomials 

Li(£) = l, L 2 (£) = £, L&)=\(. 3? 2 -l), L 4 f) = i(5^ 3 -3?). (13.3) 

Here <q through c 4 have dimension of length. Functions (13.3) and their first two £-derivatives are plotted in 
Figure 13.1. Unlike the shape functions (12.12), the L, have a clear physical meaning: 

L 1 Translational rigid body mode. 

L 2 Rotational rigid body mode. 

L 3 Constant-curvature deformation mode, symmetric with respect to £ = 0. 

L 4 Linear-curvature deformation mode, antisymmetric with respect to £ = 0. 

These properties are also shared by the standard polynomial <q + c 2 £ + of 2 + c 4 § 3 . What distinguishes the 
set (13.3) are the orthogonality properties 


Q n = / L t Ldx ~ l diag [ 1 1/3 1/5 1/7], Q 2 = / (L") r L" dx = — diag [0 0 3 25] 


Q 3 = f (L'") r L '"dx = —^diag[0 0 0 1], in which (.)'=—. 

Jo P dx 

(13.4) 

Q„ is called the covariance matrix for the n th derivative of the Legendre polynomial interpolation. The 
first-derivative covariance Qj = (L ') T L dx is not diagonal, but this matrix is not used here. 


1 They do not reach, however, the capstone level of Advanced Finite Element Methods [?]. 
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FIGURE 13.1. The Legendre polynomials and their first two?-derivatives shown over ?g[— 1, 1], Those 
interpretable as beam rigid body modes (L i and Li) in black; deformational modes (L 3 and L 4 ) in color. 


Remark 13.1. The notation (13.2)—(13.3) is FEM oriented. L\ through L 4 are called Pq through P 3 in the mathematical 
literature; e.g. Chapter 22 of the handbook [2]. The general definition for n = 0, 1 ... is 


G+i(?) ee P n (f) = 

k =0 






k =0 


(”) (?-l)" _ *(?+l)* = 



where is the binomial coefficient. Legendre polynomials are normalized by P„(l) = 1, P n (— 1) = 
also be indirectly defined by generating functions such as 


2n-2k^ n _ 2k 

(13.5) 
(—1)". They can 



]^Pn($)z" = 


k =0 


-y/1 — 2 + Z 2 


■ alternatively (?)z"=e«/o {Zyf 

k =0 


l-? 2 ). 


(13.6) 


They can also be defined through a 3-term recurrence relation: (n + 2)P„ +2 (§) — (2 n + 3) P n + i(f) + (n + 1) P„(?) = 0 
started with Po(?) = 1 and P\ (?) = ?• One important application of these polynomials in numerical analysis is the 
construction of one-dimensional Gauss integration rules: the abscissas of the n-point rule are the zeros of L n +i (?) = (f). 


§13.2.2. Generalized Stiffnesses 

The beam stiffness matrix expressed in terms of the c, is called a generalized stiffness. Denote the beam 
bending and shear rigidities by R B and R s , respectively. Then K, = K r/; + K CiS , where K cB comes from the 
bending energy and K rS from the shear energy. For the latter is its assumed that the mean shear distortion y at 
a cross section is y = T l 2 v'", where T is a dimensionless coefficient that depends on the mean-shear model 
used. Then 

R b (L") r L" dx, K cS = R s T 2 i A (L '") 71 L'” dx. (13.7) 

Jo 

In the case of a Bernoulli-Euler (BE) beam, the shear contribution is dropped: K c = K r/j . Furthermore if 
the element is prismatic, R B = El is constant. If so K r , 5 = R B and = El Q 2 , where Q 2 is the second 
diagonal matrix in (13.4). With view to future use it is convenient to differentiate between symmetric and 
antisymmetric bending rigidities R Bs and R Ba , which are associated with the responses to modes L 3 and L 4 , 
respectively. Assuming R Bs and R Ba to be uniform along the element we get 

144 R Bs 1200 R Ba 

K c = K cBs +K fBfl , K cBs = — ^diag[0 0 10], K cBa = ^ diag[0 0 0 1], (13.8) 



If shear flexibility is accounted for, the contribution K (;5 of (13.7) is kept. Assuming R s to be constant over 
the element, K r is split into 3 contributions (two bending and one shear): 


K c — + K, 


V:.S’ > 


with K fS = R s T 2 1 4 Q 3 = 


14400/?sT 2 


diag [ 0 0 0 1], 


(13.9) 
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§13.2.3. Transforming to Physical Freedoms: BE Model 

For a BE beam model, the generalized coordinates c, of (13.2) can be connected to the physical DOFs by 


V\ 
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(13.10) 


In compact form: u = G B c and c — H /; u\ with H /; = G fi '. Here 9\ = v[ and 9 2 = v 2 , which reflects the 
fundamental “plane sections remain plane” kinematic assumption of the BE model. The physical stiffness is 


K e 


H£ (K cBs + K cBa ) H b = 



12 Ra 

6 R a l 

— 12 R a 

6 R a l 

1 

6R a l 

(3 R a + R s )l 2 

-6 R a l 

(3 Ra - Rs)l 2 

£3 

-12 R a 

-6 R a l 

12 R a 

-6 R a l 


6 R a l 

(3 R a ~ Rs)l 2 

-6 R a l 

(3 R a + R s )l 2 


(13.11) 


If R s = R a = El the well known stiffness matrix (12.20) is recovered, as can be expected. The additional 
freedom conferred by (13.11) is exhibited later in two unconventional applications. 

§13.2.4. Transforming to Physical Freedoms: Shear-Flexible Model 

A shear flexible beam has mean shear distortion y = T cri 2 v'". If T is constant and u(£) interpolated by 
(13.2), v'" = 120 c 4 /l 3 . Thus y = 120T c 4 /l is constant over the element. The end rotational freedoms 
become 9\ = v[ + y and 9 2 = v' 2 + y■ Using $ = 12T to simplify the algebra, the transformations (13.10) 
change to 
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(13.12) 


In compact form, n = G,y c and c = G s u = H v uL Transforming K, of (13.9) to physical freedoms yields 
the stiffness used to construct the Timoshenko beam element in §13.2.6: 


K e 


h£k c h s = 


Rbs 

T 


"0 

0 

0 

0 " 


" 4 

21 

-4 

21 ~ 

0 

1 

0 

-1 

12 R Ba + R s <5>H 2 

21 

e 

-21 

l 1 

0 

0 

0 

0 

4£ 3 (1 + $) 2 

-4 

-21 

4 

-21 

0 

-1 

0 

1 


_ 21 

l 2 

-21 



(13.13) 


§13.2.5. Hinged Plane Beam Element 

The two-node prismatic plane BE beam element depicted in Figure 13.2 has a mechanical hinge at midspan 
(§ = 0). The cross sections on both sides of the hinge can rotate respect to each other. The top figure also 
sketches a fabrication method sometimes used in short-span pedestrian bridges. Gaps on either side of the 
hinged section cuts are filled with a bituminous material that permits slow relative rotations. 

Both the curvature k and the bending moment M must vanish at midspan. But in a element built via cubic 
interpolation of v(x), k = v” must vary linearly in both £ and x. 

Consequently the mean curvature, which is controlled by the Legendre function L 2 (shown in blue on Fig¬ 
ure 13.1) must be zero. The kinematic constraint of zero mean curvature is enforced by setting the symmetric 
bending rigidity R Bs =0 whereas the antisymmetric bending rigidity is the normal one: R Ba = El. 
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Figure 13.2. Beam element with midspan hinge. Left figure sketches a hinge fabrication method. 


Plugging into (13.11) yields 
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(13.14) 


This matrix has rank one, as it can be expected from the last (dyadic) expression in (13.14). K' has one 
nonzero eigenvalue: 6EI (4 + 1 2 )/l 3 , and three zero eigenvalues. The eigenvector associated with the nonzero 
eigenvalue pertains 2 Matrix (13.14) can be derived by more sophisticated methods (e.g., mixed variational 
principles) but the present technique is the most expedient one. 


Example 13.1. This example deals with the prismatic continuous beam shown in Figure 13.3(a). This has two spans with 
lengths ah and L, respectively, where a is a design parameter, and is subjected to uniform load q$. The beam is free, simple 
supported and fixed at nodes 1, 2 and 3, respectively. There is a hinge at the center of the 2-3 span. Of interest is the design 
question: for which a > 0 is the tip deflection at end 1 zero? 

The beam is discretized with two elements: (1) and (2), going from 1 to 2 and 2 to 3, respectively, as shown in the figure. 
The stiffnesses for elements (1) and (2) are those of (12.20) and (13.14), respectively, whereas (12.21) is used to build the 
consistent node forces for both elements. 


Assembling and applying the support conditions u 2 = D = #3 = 0, provides the reduced stiffness equations 


- 12 /a 3 

6 L/a 2 

6L/a 2 
AL 2 /a 

6L/a 2 

2 L 2 /a 

~v\ “ 

0 \ 

_ qo l 

a 

La 2 /6 

-6L/a 2 

2 L 2 Iol 

L 2 ( 4 + 3a)/a _ 

-O 2 - 

z 

- 1 

<N 

S 

1 


(13.15) 


Solving for the node displacements gives v\ = qoL 4 a(12a 2 + 9a 3 — 2)/(72£7), = qoL 3 (l — 6 a 2 — 6a 3 )/(36EI) 

and 62 = qoL 3 {l — 6a 2 )/(36EI). The equation iq = 0 is quartic in a and has four roots, which to 8 places are 
0:1 = —1.17 137898, ai = —0.52399776, 013 = 0, and 0:4 = 0.3620434. Since the latter is the only positive root, the 
solution is a = 0.362. The deflection profile for this value is pictured in Figure 13.3(b). 



(b) 



FIGURE 13.3. Beam problem for Example 13.1. (a): beam problem, (b) deflection for a = 0.362. 


2 Compare the vector in the last expression in (13.14) to the last row of Hg in (13.10) to the antisymmetric deformational 
mode L 3 . 
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§13.2 GENERALIZED INTERPOLATION 



(C) y 



A positive transverse 
shear force V produces 
a CCW rotation (+y) of 
the beam cross section 


FIGURE 13.4. Two-node Timoshenko plane beam element: (a) kinematics (when developed with cubic shape 
functions, y\ = Y 2 = Y)\ (b) M and V sign conventions; (c) concurrence of sign conventions for V and y. 


§13.2.6. Timoshenko Plane Beam Element 

As observed in §12.2.2, the Timoshenko beam model [809] includes a first order correction for transverse 
shear flexibility. The key kinematic assumption changes to “plane sections remain plane but not necessarily 
normal to the deformed neutral surface.” This is illustrated in Figure 13.4(a) for a 2-node plane beam element. 
The cross section rotation 6 differs from t/ by y. ignoring axial forces, the displacement field is analogous to 
that of the Bemoulli-Euler model (12.2) but with a shear correction: 

dv , V 

u(x, y) = — yd, v(x, y ) = v(x), with 9 =-1 - y = v + y, y = -. (13.16) 

dx GA S 

Here V is the transverse shear force, y the “shear rotation” (positive CCW) averaged over the cross section, G 
the shear modulus and A s the effective shear area. 3 The product R s = GA S is the shear rigidity. To correlate 

with the notation of §13.2.4, note that V = El v"', y = T£ 2 v'” = V/(GA S ), so T = El/{GA S £ 2 ) and 


<I> = 12T 


\2EI 

GA S i 2 


(13.17) 


This dimensionless ratio characterizes the “shear slenderness” of the beam element. 4 it is not an intrinsic 
beam property because it involves the element length. As —> 0 the Timoshenko model reduces to the BE 
model. Replacing R Bs = R Ba = El and R s = GA S = 12£7/(<M 2 ) into (13.13) yields the Timoshenko 
beam stiffness 


K e 


El 

£ 3 (1 + «&) 


12 
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-12 

61 

61 

f{A + <D) 

-61 

t 2 {2 - d>) 

-12 
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12 
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(13.18) 


If <b = 0 this reduces to (12.20). The Mathematica module TimoshenkoBeamStif f ness [Le, El, O], listed 
in Figure 13.5, implements (13.18). 


3 A concept defined in Mechanics of Materials; see e.g. Chapter 10 of Popov [656] or Chapter 12 of Timoshenko and 
Goodier [811]. A s is calculated by equating the internal shear energy j Vy = ^ V 2 /(GA S ) to that produced by the shear 
stress distribution over the cross section. For a thin rectangular cross section and zero Poisson’s ratio, A s = 5A/6. 

4 Note that in (13.8)—(13.9), 1200 R Bs /l 3 + 14400 R B T 2 /l = 1200 El (1 + 4>)/f 3 , giving a simple interpretation for <t>. 
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TimoshenkoBeamStiffness [Le_, EI_, <t>_] : =Module [ {Ke}, 

Ke=EI/(Le* (l+T) )*{{ 12/Le A 2, 6/Le,-12/Le A 2, 6/Le }, 

{ 6/Le , 4+<h, -6/Le , 2-<h }, 

{-12/Le A 2, -6/Le, 12/Le A 2,-6/Le}, 

{ 6/Le , 2-<h, -6/Le , 4+<h }}; 

Return[Ke]]; 


FIGURE 13.5. Module to produce stiffness matrix for Timoshenko beam element. 


§13.2.7. Shear and Curvature Recovery 

When using the Timoshenko beam model, the following arises during postprocessing. Suppose that the 
element node displacement vector u = u' r is given following the solution process. Recover the mean shear 
distorsion y e and the curvature k'' over the element on the way to internal forces and stresses. The problem 
is not trivial because y e is part of the rotational freedoms. The recovery process can be effectively done 
by passing first to generalized coordinates: c e = H s iT , and then to Bernoulli-Euler node displacements: 
u e BE = c e = G s Hy iT = T bt u e T , in which 
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Subtracting: u y = = (I — G B H 5 )u^ = [0 y e 0 y e ] T , Explicitly 

Y e = ^(v e i-v e 2 ) + kW + 80 


T BT — — Gfl Hy — I — 
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0 

j_ 

L e 

0 

J_ 

L e 


(13.19) 


(13.20) 


The curvature is obtained from the Bernoulli-Euler vector: k = Wn e BE , where the curvature displacement 
matrix B' is that given in the previous Chapter. 

§13.2.8. Generalized and Physical Loads 

Assume that the Legendre-modeled beam is loaded by a distributed lateral force specified as a polynomial 


<?(£) — a Q + a l t + a 2 £ +<T? 3 - 


(13.21) 


This is converted to generalized consistent loads through the integral 


_ 1 

~L i(£)~ 


CIq -f- 0.2/33 

5 = 5* j 

l 2 (?) 

l 3 (?) 

d%=l 

0\/3 o 3 /5 
2a 2 /l5 

l 4 (?) 


2a 2 /35 


(13.22) 


For the BE model, transforming (13.22) to physical coordinates via f = W T B f„ with the H fi of (13.10) yields 
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The first term reproduces the uniform load result (12.21) of the previous Chapter if q = a 0 , others zero. 
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For the Timoshenko (shear flexible) beam, using the H s of (13.12) we get 
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(13.24) 

The terms pertaining to ao and a 2 , which are centrosymmetric, do not depend on <t> and agree with those for 
the BE beam model. 

Often q is specified in terms of node values, e.g., a linearly varying load q = q\{\ — §)/2 + q 2 (1 + §)/2. 
This can be mapped to (13.21) by setting a Q = \{q\ + q 2 ), + = §(— q\ + q 2 ), a 2 = a 2 = 0. 



Figure 13.6. Beam for Example 13.2: cantilever discretized with two Timoshenko beam elements. 


Example 13.2. Consider the prismatic cantilever beam of length L pictured in Figure 13.6(a). It is subject to two point 
loads as shown. Shear flexibility is to be accounted for using the Timoshenko model. The bending and shear rigidities E I zz 
and GA S are constant along the span. The objective is to find deflections, curvatures and shear distortions and associated 
bending moments and shear forces. 

It is sufficient to discretize the beam with two Timoshenko beam elements oflength L /2 as shown in the figure. The stiffness 
matrices for both elements are given by (13.18), in which L e = L/2and 4> = 12 E l zz /(G A s (L/2) 2 ) = 48 E I ZZ /(GA S L 2 ). 
The master stiffness equations are 
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Setting the displacement B.C. tq = 9\ = 0 and solving yields 

U = I37^[ 0 0 i 1 fa 122 + 0) ■ 

The mean element shear distortions are calculated from (13.26) using (13.20). This gives 


y (1) =0, 


48 E I zz G A s 


(13.25) 


(13.26) 


(13.27) 
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Figure 13.7. Beam supported by a bed of springs. Continuification of this configuration 
leads to the Winkler foundation model treated in this subsection. 


The element-level Bernoulli-Euler node displacements are obtained from (13.26) on subtracting the shear distortions (13.27) 
from the rotations: 
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(13.28) 


Note that 9 BZ 
recovered as 


(2) | 

1 fi = 1 + j 2 3>, the “kink” being due to the shear distortion jump at node 2. The curvatures are now 


4D 


B ( 1 ) u^ = 
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B (2) u^' = 


PL 
4 EL, 


-d-? <2> ) 


(13.29) 


The transverse shear force resultant and bending moment are easily recovered as V‘ = G A s y (e) and Mt = E I zz K e , 
respectively. The results are drawn in Figure 13.6(b,c). 


§13.2.9. Beam on Elastic Supports 


Sometimes beams, as well as other structural members, may be supported elastically along their span. Two 
common configurations that occur in structural engineering are: 

(i) Beam resting on a continuum medium such as soil. This is the case in foundations. 

(ii) Beam supported by discrete but closely spaced flexible supports, as in the “bed of springs” pictured in 
Figure 13.7. This occurs in railbeds (structurally rails are beams supported by crossties) and some types 
of grillworks. 

The Winkler foundation is a simplified elastic-support model. It is an approximation for (i) because it ignores 
multidimensional elasticity effects as well as friction. It is a simplification of (ii) because the discrete nature 
of supports is smeared out. It is nonetheless popular, particularly in foundation and railway engineering, when 
the presence of physical uncertainties would not justify a more complicated model. Such uncertainties are 
inherent in soil mechanics. 

The Winkler model may be viewed as a “continuification” of case (ii). Take a beam slice going from x to 
x + dx. The spring-reaction force acting on the beam is taken to be df F = —k F v(x) dx. Here v(x) is the 
transverse deflection and k F the Winkler foundation stiffness, which has dimension of force per length-squared. 
Force df F has the opposite sign of v(x), pushing up if the beam moves down and pulling down if it moves up. 
Beam-foundation separation effects that may occur in case (i) are ignored here because that would lead to a 
nonlinear contact problem. 

The internal energy stored in the dx slice of Winkler springs is — 1/2 v df F = 1/2 k r v 2 dx. Consequently the 
effect of elastic supports is to modify the internal energy U B of the beam element so that it becomes 


U e = U e B + U e F , with U e F = \ f k F v 2 dx. (13.30) 

Jo 
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§13.3 INTERPOLATION WITH HOMOGENEOUS ODE SOLUTIONS 


Therefore the total stiffness of the element is computed by adding the foundation stiffness to the beam stiffiiess. 
Care must be taken, however, that the same set of nodal freedoms is used. This is best handled by doing the 
generalized stiffness first, and then using the appropriate generalized-to-physical transformation. If the 
transverse deflection v is interpolated with (13.2) as v — L c, the generalized Winkler foundation stiffness for 
constant k F is 

K cF =k F f L r Ldx = k F Q 0 , (13.31) 

Jo 


where Q 0 is the first diagonal matrix in (13.4). This holds regardless of beam model. Now if the member resting 
on the foundation is modeled as a BE beam, one picks H B of (13.10) as generalized-to-physical transformation 
matrix to get 


K 


e 

F 



156 

22 i 

54 

-13 1 

k F i 

22£ 

41 2 

13£ 

—3£ 2 

420 

54 

\3£ 

156 

—22£ 


— 13 £ 

-3£ 2 

—221 

4£ 2 


(13.32) 


If instead the supported member is modeled as a Timoshenko beam, one picks H v of (13.12) to get 


K e F = k F H T s Q 0 H s 

" 4(78+147cb+70<& 2 ) £ (44+774>+35<I> 2 ) 4(27+634>+354> 2 ) -£ (26+63 4>+354> 2 ) " 
k F l f (44+774>+35<T 2 ) f 2 (8+14<fi+7<i> 2 ) l (26+63<J>+354> 2 ) -f 2 (6+144>+74> 2 ) 

~ 840(l+<f>) 2 4(27+634>+354> 2 ) £ (26+63$+354> 2 ) 4(78+147<b+70cfi 2 ) -l (44+77<3>+350 2 ) 

—£(26+630+354> 2 ) —€ 2 (6+140+7 <t* 2 ) —f (44+774>+35<T 2 ) € 2 (8+14<t>+7<T* 2 ) 

(13.33) 

The module TimoshenkoWinklerStif fness [Le ,kF, 4>] listed in Figure 13.8 implements the stiffness 
(13.33). To get the BE-beam Winkler stiffness (13.32), invoke with = 0. Examples of use of this module 
are provided in §13.3.1. 


TimoshenkoWinklerStiffness[Le 

kF_, 4>_] : =Module [ {KeW}, 

KeW={ {4* (78+147*4*70*4> A 2) , 

Le* (44+77*4*35*4^2), 

4* (27 + 63*4*35*4> A 2) , 

-Le*(26+63*4*35*4^2)}, 

{Le*(44+77*4*35*4^2), 

Le A 2*(8+14*4*7*4^2), 

Le*(26+63*4*35*4^2), 

-Le A 2* (6+14*4*7*4^2)}, 

{4*(27+63*4*35*4^2), 

Le* (26+63*4*35*4^2), 

4*(78+147*4*70*4^2), 

-Le*(44+77*4*35*4^2)}, 

{-Le*(26+63*4*35*4^2), 

-Le A 2*(6+14*4*7*4^2), 

-Le* (44+77*4*35*4^2), 

Le A 2*(8+14*4*7*4^2)}}* 

kF*Le/ (840* (l+4>) A 2) ; 

Return[KeW]]; 


FIGURE 13.8. Stiffness matrix module for a Winkler foundation supporting a Timoshenko beam element. 


§13.3. Interpolation with Homogeneous ODE Solutions 

For both BE and Timoshenko beam models, the Legendre polynomials Lff) through L 4 (f ) are exact solutions 
of the homogeneous, prismatic, plane beam equilibrium equation El d A v/dx 4 = 0. When used as shape 
functions in the generalized interpolation (13.2), the resulting stiffness matrix is exact if the FEM model is 
loaded at the nodes, as further discussed in §13.6. The technique can be extended to more complicated one¬ 
dimensional problems. It can be used to derive exact stiffness matrices if homogeneous solutions are available 
in closed form, and are sufficiently simple to be amenable to analytical integration. The following subsection 
illustrates the method for a BE beam resting on a Winkler elastic foundation. 
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Chapter 13: ADVANCED ONE-DIMENSIONAL ELEMENTS 


§13.3.1. Exact Winkler/BE-Beam Stiffness 


Consider again a prismatic, plane BE beam element resting on a Winkler foundation of stiffness k F , as pictured 
in Figure 13.7. The governing equilibrium equation for constant El > 0 and kp > 0 is El d 4 v/dx 4 + k F v — 
q(x). The general homogeneous solution over an element of length l going from x = 0 to x = i is 


v(x) = e i; (c i sin£ +C2 cosf) +e ( (c 3 sin£ + C 4 cost;), with £ = X x /Z and x 



(13.34) 


Here the c, are four integration constants to be determined from four end conditions: the nodal degrees of 
freedom v\, v\, v 2 and v' 2 . These constants are treated as generalized coordinates and as before collected 
into vector c = [ c \ c 2 c 3 c 4 ] 7 . The solution (13.34) is used as generalized interpolation with <?'■ sin £ through 
e~^ cos £ as the four shape functions. Differentiating twice gives v' = dv/dx and v" = d 2 v/dx 2 . The TPE 
functional of the element in terms of the generalized coordinates can be expressed as 


n e c = [ (\EIW’?+ \k F v 2 -q 0 v) dx= ±c T (K cB +K cF )c-c r f c . (13.35) 

Jo 


This defines and K cf as generalized element stiffnesses due to beam bending and foundation springs, 
respectively, whereas f, is the generalized force associated with a transverse load q(x). The nodal freedoms 
are linked to generalized coordinates by 


V\ 


Vl 

0i 


0 1 0 

X/Z X/Z _X/Z 

e x sin x e x cos x e~ x sin x 

Xe x (cos x+sin x) x^fcos X^ sin x) x g_x (cos x— sinx) — X e 

t l t 


1 

-X/Z 
e x cos x 
x (cos x + sin x) 
T 


Cl 

Cl 

C 3 


_C 4 J 
(13.36) 


In compact form this is iT” = G F c. Inverting gives c = u' with H/r = G F '. The physical stiffness is 
K' = K) ; + K' /r with = H 7 K, /; H F and Kg = W] K cF Hg. The consistent force vector is P = H 7 - f r . 
Computation with transcendental functions by hand is unwieldy and error-prone, and at this point it is better 
to leave that task to a CAS. The Mathematica script listed in Figure 13.9 is designed to produce Kg, Kg and 
f for constant El and k F , and uniform transverse load q(x) = q Q . The script gives 



B 1 B 2 

b 5 

-b 4 ~ 


F 1 F 2 

f 5 

-f 4 ~ 


“ /l ■ 

El x 

Bi 

b 4 

Be 

K - kFi 

f 3 

f 4 

f 6 

ye doZ 

fl 

At 3 g 1 


Bi 

-b 2 

’ F 16x 3 g 2 


F\ 

-f 2 

’ ^ 2 
X^-g 

fl 


symm 


b 3 


symm 


f 3 


-fl. 


(13.37) 
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§13.3 INTERPOLATION WITH HOMOGENEOUS ODE SOLUTIONS 


ClearAll[El,kF,Le,%,qO,x]; g=2-Cos[2 *%]-Cosh[2*%]; 

Nf={Exp[ X*x/Le]*Sin[X*x/Le],Exp[ %*x/Le]*Cosft*x/Le], 

Exp[-x*x/Le]*Sinft*x/Le],Exp[~X*x/Le]*Cosft*x/Le]}; 

Nfd=D[Nf,x]; Nfdd=D[Nfd,x]; 

KgF=kF*Integrate[Transpose[{Nf)].{Nf},{x,0,Le}]; 
KgB=EI*Integrate[Transpose[{Nfdd}].{Nfdd},{x,0,Le}]; 
fg= qO*Integrate[Transpose[{Nf}],{x,0,Le}]; 

{KgF,KgB,fg}=Simplify[{KgF,KgB,fg}]; 

Print["KgF=",KgF//MatrixForm]; Print["KgB=",KgB//MatrixForm]; 
GF=Simplify[{Nf/.x->0,Nfd/.x->0,Nf/.x->Le,Nfd/.x->Le}]; 
HF=Simplify[Inverse[GF]]; HFT=Transpose[HF]; 

Print["GF=",GF//MatrixForm]; Print["HF=",HF//MatrixForm]; 
facB=(EI*x/Le A 3)/(4*g A 2); facF=(kF*Le)/(16*X A 3*g A 2) ; 
KeB=Simplify[HFT.KgB.HF]; KeBfac=FullSimplify[KeB/facB]; 

Print["KeB=",facB," * ",KeBfac//MatrixForm] 

KeF=Simplify[HFT.KgF.HF]; KeFfac=FullSimplify[KeF/facF]; 

Print["KeF=",facF," * ",KeFfac//MatrixForm]; 

facf=(qO*Le)/(X A 2*g); fe=Simplify[ExpToTrig[HFT.fg]]; 

fefac=Simplify[fe/facf]; Print["fe=",facf," * ",fefac]; 


FIGURE 13.9. Script to produce the exact Winkler-BE beam stiffness matrix and consistent force vector. 


BEBeamWinklerExactStiffness[Le_,EI_,kF_,q0_]:=Module[{B1,B2,B3,B4,B5,B6, 
F1,F2,F3,F4,F5,F6,f1,f2,facB,facF,facf,KeB,KeF,fe,X}, 

X=PowerExpand[Le*((kF/(4*EI)) A (1/4))]; 

B1 =2*X a 2*(-4*Sin[2*X]+Sin[4*X]+4*Sin[x]*(Cos[X]*Cosh[2*X]+ 

8*X*Sin[X]*Sinh[X] A 2)+2*(Cos[2*X]-2)*Sinh[2*X]+Sinh[4*X]); 

B2 =2*Le*X*(4*Cos[2*X]-Cos[4*X]-4*Cosh[2*X]+ Cosh[4*X]- 
8*X*Sin[2*X]*Sinh[X] A 2+8*X*Sin[X] A 2*Sinh[2*X]); 

B3 =-(Le A 2*(8*X*Cos[2*X]-12*Sin[2*X]+Cosh[2*X]*(6*Sin[2*X]-8*X)+3*Sin[4*X]+ 

2*(6-3*Cos[2*X]+4*X*Sin[2*X])*Sinh[2*X]-3*Sinh[4*X])); 

B4 =-4*Le*X*(X*Cosh[3*X]*Sin[X]~X*Cosh[x]*(-2*Sin[X]+Sin[3*X]) + (X* (Cos [X] + 

Cos[3*x])+Cosh[2*x]*(-2*X*Cos[X]+4*Sin[X])+2*(-5*Sin[X]+Sin[3*X]))*Sinhft ]) ; 
B5 =-4*X A 2* (2*Cos[X]*(-2+Cos[2*X]+Cosh[2*X])*Sinh[X]+Sin[3*X]* 

(Cosh[X]-2*X*Sinh[X])+Sin[X]*(-4*Cosh[X]+Cosh[3*X]+2*X*Sinh[3*X])); 

B6 =2*Le A 2*(Cosh[3*X]*(-2*X*Cos[X]+3*Sin[X])+Cosh[x]*(2*X*Cos[3*X]+3*(Sin[3*X]- 
4*Sin[X]))+(9*Cos[X]-3*Cos[3*X]-6*Cos[X]*Cosh[2*x]+16*X*Sin[X])*Sinh[X]); 

FI =2*X a 2*(-32*X*Sin[x] A 2*Sinh[X] A 2+6* (-2+Cos [2*X])* 

(Sin[2*X]+Sinh[2*X])+6*Cosh[2*X]*(Sin[2*X]+Sinh[2*X])); 

F2 =2*Le*X* (4*Cos[2*X]-Cos[4*Xl-4*Cosh[2*X]+Cosh[4*X] + 

8*X*Sin[2*X]*Sinh[X] A 2-8*X*Sin[X] A 2*Sinh[2*X]); 

F3 =Le A 2*(8*X*Cos[2*X]+4*Sin[2*X]-2*Cosh[2*X]*(4*X+Sin[2*X])-Sin[4*X]+ 
2*(Cos[2*X]+4*X*Sin[2*X]-2)*Sinh[2*X]+Sinh[4*X]); 

F 4 =4*Le*X*(X*Cosh[3*X]*Sin[x]~X*Cosh[X]*(-2*Sin[X]+Sin[3*X]) + (X*Cos[X] + 

X*Cos[3*X]+10*Sin[X]-2*Cosh[2*X]*(X*Cos[X]+2*Sin[X])-2*Sin[3*X])*Sinh[X]); 
F5 =-4*X A 2*(6*Cos[X]*(-2+Cos[2*X]+Cosh[2*X])*Sinh[x]+Sin[3*X]* 

(3*Cosh[x]+2*X*Sinh[x])+Sin[x]*(-12*Cosh[x]+3*Cosh[3*x]-2*x*Sinh[3*x])); 

F6 =-2* Le A 2*(-(Cosh[3*X]*(2*x*Cos[X]+Sin[x]))+Cosh[x]* (2*x*Cos [3*x] + 

4*Sin[x]-Sin[3*X])+(Cos[3*X]+Cos[X]*(2*Cosh[2*X]-3)+16*X*Sin[x])*Sinh[X]); 
f1=2*X*(Cosh[X]"Cos[X])*(Sin ft]-Sinh[X]); f2=-(Le*(Sin[x]-Sinh[X]) A 2); 
g=2-Cos[2 *x]-Cosh[2*x]; facf=(qO*Le)/(x A 2*g); 

facB=(EI*x/Le A 3)/(4*g A 2); facF=(kF*Le)/(16*X A 3*g A 2); 

KeB=facB*{{B1,B2,B5,-B4},{ B2,B3,B4,B6},{B5,B4,Bl,-B2),{-B4,B6,-B2,B3}}; 
KeF=facF*{{FI,F2,F5,-F4},{ F2,F3,F4,F6},{F5,F4,FI,—F2},{-F4,F6,-F2,F3}}; 
fe=facf*{f1,f2,f1,-f2}; Return[{KeB,KeF,fe}]]; 


FIGURE 13.10. Module to get the exact BE-Winkler stiffness and consistent load vector. 


13-13 
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y, vA 


El constant 


IP : Load case (I) 

▼ q Q : Load case (II) 



FIGURE 13.11. Example: fixed-fixed beam on Winkler elastic foundation. 


in which 

g = 2 — cos 2/ — cosh 2x, 

B | = 2x 2 (—4 sin2x+ sin4x+4 sin x (cos x cosh 2x+8x sin x sinlr x)+2(cos 2x—2) sinh2x+ sinh4x), 

Bi = 2^x(4cos2x— cos4x—4cosh2x + cosh4x—8 x sin2x sinh 2 x+8x sin 2 X sinh2x), 

B 3 = —(£ 2 (8x cos2x —12sin2x+cosh2x(6sin2x—8x)+3sin4x + 

2(6—3cos2x+4x sin2x) sinh2x—3 sinh4x)), 

B 4 = —4£x(x cosh3x sinx—X cosh x(—2sinx+ sin3x)+(x( cos X + cos 3x) 

+ cosh2x(—2x cos x+4sin x)+2(—5 sinx+ sin3x)) sinhx), 

B 5 = —4x 2 (2cos x (—2+ cos 2x+ cosh 2x) sinh x + sin 3x (cosh x— 2x sinh x) 

+ sinx(—4coshx+cosh3x+2x sinh3x)) 

B 6 = 2^ 2 (cosh3x(—2x cos x+3 sin x)+ cosh x(2x cos 3x+3(—4 sin x+ sin 3x)) 

+(9cosx—3cos3x—6cosx cosh2x+16x sinx) sinhx) 

F\ = 2x 2 (—32x sin 2 x sinh 2 x+6(—2+cos2x)(sin2x+ sinh2x)+6cosh2x(sin2x + sinh2x)), 

F 2 = 2^x(4cos2x— cos4x—4cosh2x + cosh4x+8x sin2x sinh 2 x~ 8x sin 2 x sinh2x), 

F 3 = t! 2 (8x cos2x+4sin2x—2cosh2x(4x + sin2x)— sin4x+2(cos2x+4x sin2x—2) sinh2x+ sinh4x), 
F a = 4 lx(x cosh3x sinx—X cosh x(sin3x~2sinx)+(x cos X+X cos3x+10sinx 
—2cosh2x(x cos x+2 sin x)—2sin3x) sinh x), 

F 5 = —4x 2 (6cosx(—2+cos2x+cosh2x)sinhx + sin3x(3coshx+2x sinhx), 

+ sin x(—12 cosh x+3 cosh 3x—2x sinh3x)) 

F 6 = — 2l 2 {— (cosh3x(2x cosx +sinx))+coshx(2x cos3x+4sinx — sin3x) 

+(cos 3x + cos x (2 cosh 2x — 3)+16x sin x) sinh x), 

/i = 2x(coshx - cosx)(sinx - sinhx), fi = -^(sinx - sinhx) 2 . 

(13.38) 

These expressions are used to code module BEBeamWinklerExactStiffness[Le,EI,kF,qO], which is 
listed in Figure 13.10. 

Example 13.3. A fixed-fixed BE beam rests on a Winkler foundation as shown in Figure 13.11. The beam has span 2L, and 
constant El. The Winkler foundation coefficient kp is constant. As usual in foundation engineering we set 

k F = EIk A /L 4 , (13.39) 
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§13.3 INTERPOLATION WITH HOMOGENEOUS ODE SOLUTIONS 


Table 13.1 - Results for Example of Figure 13.11 at Selected 1 Values 



Load case (I): Central Point Load 

Load case (II): Line Load Over Right Half 

X 

exact C i 

N e = 2 

II 

7 

N e = 8 

exact Cjj 

II 

=5 

II 

7 

N e = 8 

0.1 

0.999997 

0.999997 

0.999997 

0.999997 

0.999997 

0.999997 

0.999997 

0.999997 

1 

0.970005 

0.969977 

0.970003 

0.970005 

0.968661 

0.969977 

0.968742 

0.968666 

2 

0.672185 

0.668790 

0.671893 

0.672167 

0.657708 

0.668790 

0.658316 

0.657746 

5 

0.067652 

0.049152 

0.065315 

0.067483 

0.041321 

0.049152 

0.041254 

0.041317 

10 

0.008485 

0.003220 

0.006648 

0.008191 

0.002394 

0.003220 

0.002393 

0.002395 

100 

8.48x 10~ 6 

3.23 x 10~ 7 

8.03 x 10~ 7 

1.63x 10~ 6 

2.40x 10~ 7 

3.23 x 10~ 7 

2.62x 10~ 7 

2.42 x 10" 7 


where X is a dimensionless rigidity to be kept as parameter. 5 The beam is subjected to two load cases: (I) a central point 
load P at x = L , and (II) a uniform line load qo over the right half .* > L. See Figure 13.11. 

All quantities are kept symbolic. The focus of interest is the deflection vc at midspan C (x = L). For convenience this 
is rendered dimensionless by taking v ! c (X) = CiiXjv^ and vjJ (X) = Cu(X)v I ( J 0 for load cases (I) and (II)), respectively. 
Here = — P L? / (24EI) and v{J Q = —qoL 4 /(48EI) are the midspan deflections of cases (I) and (II) for X = 0, that is, 
kf = 0 (no foundation). The exact deflection factors for this model are 


Ci(X) = 


Cn(X) = 


6\/2 cos ~JlX + cosh \/2X — 2 13 4 137 8 

k 3 sin \[7.X + sinh \f7.X 420 138600 

48 (cosk/\/2— cosh k/\/2) (sin k/\/2— sinhk/-^) i 163 , 4 i 20641 i8 

F sin V2k +sinh72k = 1 “ 5040^ + 19958400^ + 


(13.40) 


Both load cases were symbolically solved with two exact elements of length L produced by the module of Figure 13.10. 
As can be expected, the answers reproduce the exact solutions (13.40). Using any number of those elements would match 
(13.40) as long as the midspan section C is at a node. Then both cases were solved with 2,4, and 8 elements with the stiffness 
(13.32) produced by cubic polynomials. The results are shown in a log-log plot in Figure 13.12. Results for selected values 
of X are presented in Table 13.1. 

As can be seen, for a “soft" foundation characterized by X < 1, the cubic-polynomial elements gave satisfactory results and 
converged quickly to the exact answers, especially in load case (II). As X grows over one, the deflections become rapidly 
smaller, and the polynomial FEM results exhibit higher relative errors. On the other hand, the absolute errors remain small. 
The conclusion is that exact elements are only worthwhile in highly rigid foundations (say X > 5) and then only if results 
with small relative error are of interest. 



Figure 13.12. Log-log plots of C/(X) and C//(k) for Example of Figure 13.11 over range X e [0.1, 100]. 


5 Note that X is a true physical parameter, whereas x is discretization dependant because it involves the element length. 
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Remark 13.2. To correlate the exact stiffness and consistent forces with those obtained with polynomial shape functions 
it is illuminating to expand (13.37) as power series in /. The rationale is that as the element size t, gets smaller, / = 
l \J~k~p/(AEI) goes to zero for fixed El and k F . Mathematica gives the expansions 

K e B =K B0 + x 8 K B8 + x l2 K B12 + ..., K e r =K F0 + x 4 K F4 + x 8 K F8 + ..., f = f 0 + X 4 U + ■ ■ ■, (13.41) 


in which 


K^ 0 = eqn (12.20) , 


K' 


K 


El 


B 12 


5959453500 £ 3 


~ 25488 53521! 23022 -5043f 

El 5352£ 1136€ 2 5043£ -1097£ 2 

43659001 3 23022 5043£ 25488 -5352f 

_ —50431! — 10971! 2 -5352 i 1136£ 2 
113504€ 522090 -11263 If 

243841! 2 11263 If -24273f 2 

112631t 528960 -113504t 


' 528960 
113504t 
522090 

— 11263 If —24273f 2 -113504t 24384f 2 


K e F0 = eqn (13.32), 


rre 

F4 “ 2EI B8 ’ 


Ks = - 


3 k F l 4 




«S=^[6f 6 -t 


qpl 

4 5040 


[14 3f 14 —3f] y 


SEI Bl2 ' 0 12 

(13.42) 

Thus as x Owe recover the stiffness matrices and force vector derived with polynomial shape functions, as can be 
expected. Note that Kgo and Kfo decouple, which allows them to be coded as separate modules. On the other hand the 
exact stiffnesses are coupled if x > 0. The foregoing expansions indicate that exactness makes little difference if x < 1 • 


§13.4. Equilibrium Theorems 


One way to get high performance mechanical elements is to use equilibrium conditions whenever possible. 
These lead to flexibility methods. Taking advantage of equilibrium is fairly easy in one space dimension. It is 
more difficult in two and three, because it requires advanced variational methods that are beyond the scope of 
this book. This section surveys theorems that provide the theoretical basis for flexibility methods. These are 
applied to ID element construction in §13.5. 

§13.4.1. Self-Equilibrated Force System 


First we establish a useful theorem that links displacement and force transformations. Consider a FEM dis¬ 
cretized body such as that pictured in Figure 13.13(a). The generic potato intends to symbolize any discretized 
material body: an element, an element assembly or a complete structure. Partition its degrees of freedom into 
two types: r and 5. The s freedoms (5 stands for suppressed or supported) are associated with a minimal set of 
supports that control rigid body motions or RBMs. The r freedoms (r is for released) collect the rest. In the 
figure those freedoms are shown collected at invididual points P s and P, for visualization convenience. Node 
forces, displacements and virtual displacements associated with those freedoms are partitioned accordingly. 
Thus 



u = 



<5u = 


' <5u s ' 
. <5u,.. 


(13.43) 


The dimension n s of f v , u v and <5u v is 1, 3 and 6 in one-, two- and three-dimensional space, respectively. 
Figure 13.13(b) shows the force system {f s , f,.} undergoing virtual displacements, which are exaggerated for 
visibility. 6 Consider now the “rigid + deformational” displacement decomposition u = Gu v + d, in which 
matrix G (of appropriate order) represents a rigid motion and d are deformational displacements. Evaluating 
this at the r freedoms gives 

u,. = G,. u s + d, , Su,. = G,. <5u v + <5d,., (13.44) 


The first decomposition in (13.44), being linear in the actual displacements, is only valid only in geometrically 
linear analysis. That for virtual displacements is valid for a much broader class of problems. 


6 Under virtual displacements the forces are frozen for application of the Principle of Virtual Work. 
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§13.4 EQUILIBRIUM THEOREMS 



Figure 13.13. Body to illustrate equilibrium theorems. Nodal freedoms classified into supported 
( 5 ) and released (r), each lumped to a point to simplify diagram, (a) Self equilibrated node force 
system, (b) Force system of (a) undergoing virtual displacements; grossly exaggerated for visibility. 


If the supported freedom motion vanishes: u, = 0, then u, = d,. Thus d, represents a relative displacement 
of the unsupported freedoms with respect to the rigid motion Gu v , and likewise for the virtual displacements. 
Because a relative motion is necessarily associated with deformations, the alternative name deformational 
displacements is justified. 

The external virtual work is SW — <5W V + SW r = <5uJ f s + <5u ; r f, . If the force system in Figure 13.13(a) is in 
self equilibrium and the virtual displacements are imparted by rigid motions <5d r = 0 and <5u,. = G,. <5u s , the 
virtual work must vanish: <5IV = <5u 7 f v + SuJ G] f,. = <5uJ (f v + G'. f,. ) = 0. Because the <5u s are arbitrary, 
it follows that 

f, + G ( r f, = 0, f s = —Gj f r . (13.45) 

These are the overall static equilibrium equations of a discrete mechanical system in self equilibrium. Some¬ 
times it is useful to express the foregoing expressions in the complete-vector form 


u 




<5u = 


' Su s 1 

. <5u r J 



<$u s + 





Relations in (13.46) are said to be reciprocal. 1 


Remark 13.3. The freedoms in u v are virtual supports, chosen for convenience in flexibility derivations. They should not 
be confused with actual or physical supports. For instance Civil Engineering structures tend to have redundant physical 
supports, whereas aircraft or orbiting satellites have none. 

§13.4.2. Handling Applied Forces 

Consider now a generalization of the previous scenario. An externally applied load system of surface or body forces, not 
necessarily in self equilibrium, acts on the body. For example, the surface tractions pictured in Figure 13.14(a). To bring 
this under the framework of equilibrium analysis, a series of steps are required. 

First, the force system is replaced by a single resultant q, as pictured in Figure 13.14(b). 8 The point of application is P q . 
Equilibrium is restored by introducing node forces q,. and q v at the appropriate freedoms. The overall equilibrium condition 
is obtained by putting the system {q, q r , q s | through rigid-motion virtual displacements, as pictured in Figure 13.14(c). 
Point P q moves through 5u 9 , and G evaluated at P q is G ? . The virtual work is & W = <$uf (q A + G, r q,. + G 7 q) = 0 
whence 

q .5 + Gj q,. + G T q q = 0. (13.47) 

If (13.47) is sufficient to determine q v and q,., the load system of Figure 13.14(a) can be effectively replaced by the nodal 
forces —q s and —q r , as depicted in Figure 13.14(d). These are called the equivalent node forces. But in general (13.47) 
is insufficient to fully determine q s and q,. . The remaining equations to construct the equivalent forces must come from a 
theorem that accounts for the internal energy, as discussed in §13.4.3. 


7 If the model is geometrically nonlinear, the first form in (13.46) does not hold. 

8 Although the figure shows a resultant point force, in general it may include a point moment that is not shown for simplicity. 
See also Remark 13.3. 
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FIGURE 13.14. Processing non-self-equilibrated applied loads with flexibility methods, (a) Body 
under applied distributed load, (b) Substitution by resultant and self equilibration, (c) Deriving overall 
equilibrium conditions through the PVW. (d) Replacing the applied loads by equivalent nodal forces. 


Adding (13.45) and (13.47) gives the general overall equilibrium condition 

fv + q.j + G, r (f, + q,.) + G T q q = 0, (13.48) 

which is applied in §13.4.4 to the recovery of supported freedoms. 

Remark 13.4. The replacement of the applied force by a resultant is not strictly necessary, as it is always possible to write 
out the virtual work by appropriately integrating distributed effects. The resultant is primarily useful as an instructional 
tool, because matrix G q is not position dependent. 

Remark 13.5. Conditions (13.45) and (13.47), which were derived through the PVW, hold for general mechanical systems 
under mild reversibility requirements [828,§231], including geometric nonlinearities. From now on we restrict attention to 
systems linear in the actual displacements. 

§13.4.3. Flexibility Equations 

The first step in FEM equilibrium analysis is obtaining discrete flexibility equations. The stiffness equations 
introduced in Chapter 2 relate forces to displacements. At the element level they are f = K e u . By definition, 
flexibility equations relate displacements to forces: u e = F <> F, where F e is the element flexibility matrix. So 
the expectation is that the flexibility can be obtained as the inverse of the stiffness: F c = (K e )~ ! . Right? 

Wrong. Recall that K e for a disconnected free-free element is singular. Its ordinary inverse does not exist. 
Expectations go up in smoke. The same difficulty holds for a superelement or complete structure. 

To get a conventional flexibility matrix 9 it is necessary to remove all rigid body motions in advance. This can 
be done through the virtual supports introduced in §13.4.1. The support motions u s are fixed, say u s = 0. 
Flexibility equations are sought between what is left of the kinematics. Dropping the element superscript for 
brevity, for a linear problem one gets 

F rr f r = d,. (13.49) 

Note that u, does not appear: only the deformational or relative displacements. To recover u, it is necessary 
to release the supports, but if that is naively done F rr ceases to exist. This difficulty is overcome in §13.4.4. 


9 In the FEM literature it is often called simply the flexibility. The reason is that for a long time it was believed that getting 
a flexibility matrix required a supported stmcture. With the recent advent of the free-free flexibility (see Notes and 
Bibliography) it becomes necessary to introduce a “deformational” or “conventional” qualifier. 
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§13.4 EQUILIBRIUM THEOREMS 


There is another key difference with stiffness methods. The DSM assembly procedure covered in Chapter 3 
(and extended in Chapter 27 to general structures) does not translate into a similar technique for flexibility 
methods. Although it is possible to assemble flexibilities of MoM elements, the technique is neither simple 
nor elegant. And it becomes dauntingly complex when tried on continuum-based elements [251]. 

So one of the main uses of flexibility equations today is as a stepping stone on the way to derive element 
stiffness equations, starting from (13.49). The procedural steps are explained in §13.4.4. But how should 
(13.49) be derived? There are several methods but only one, based on the Total Complementary Potential 
Energy (TCPE) principle of variational mechanics is described here. 

To apply TCPE, the complementary energy II* of the body must be be expressed as a function of the nodal 
forces f,. . For fixed supports (u s = 0) and a linear system, the functional can be expressed as 

n*(f r ) = U*(f r ) - f r d, = \f r F„ f,. + f, r b, - f, r d, + n*. (13.50) 


Here U* is the internal complementary energy, also called the stress energy by many authors, e.g., [353], b, is a 
term resulting from loading actions such as as thermal effects, body or surface forces, and IT,* is independent of 
f, . Calculation of U* in ID elements involves expressing the internal forces (axial force, shear forces, bending 
moments, torque, etc.) in terms of f, from statics. Application examples are given in the next section. 10 The 
TCPE principle states that IT is stationary with respect to variations in f,. when kinematic compatibility is 
satisfied: 

an* 

—— = F„.f r + b, — d, = 0. whence d,. = F rr f,. + b,.. (13.51) 

af,. 

By hypothesis the deformational flexibility F rr is nonsingular. Solving for f, gives the deformational stiffness 
equations 

f, = K,,d r - q,., with K,. r = F“' and q r = K,,b, . (13.52) 

The matrix K,.,. is the deformational stiffness matrix, whereas q ( . is the equivalent load vector. 

§13.4.4. Rigid Motion Injection 

Suppose that F,.,. and q ( . of (13.52) have been found, for example from the TPCE principle (13.51). The goal 
is to arrive at the free-free stiffness equations, which are partitioned in accordance with (13.43) as 



(13.53) 


To justify the presence of K,.,. and q ( . here, set u s = 0, whence u r = d,. . Consequently the second equation 
reduces to f, = K, r d r — q ( . which matches (13.52). Inserting f v and f,. into (13.48) yields 


(Ky.v + G, r K ( -s)Ui + (K sr + G 1 K,.,.)u,. + [q s + G' 7 q, + G 7 q] — 0, 


(13.54) 


and replacing u r = G, u, + d, , 


[K, 


+ g; 7 k. 


+ K,G, 


+ G, r K,, 


G, ] Uj + [K. s 


+ G, r K,, 


] d, + [q v 


+ G, r q,. + G T q q] = 0. (13.55) 


Because u s , d,. and q can be arbitrarily varied, each bracket in (13.55) must vanish identically, giving 


q s = —G, 7 q, - G T q q, K w . = -G, 7 K„, K„ = Kj r = -K rr G r , 

K vv = — G, r K,. v — K V , G, — G'j. K,.,.G r = G ( 7 K, , G,. + G 7 K, , G, — G 1 K„. G r = G 7 K„. G,. 


(13.56) 


10 For 2D and 3D elements the process is more delicate and demands techniques, such as hybrid variational principles, that 
lie beyond the scope of this material. 
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Inserting these into (13.53) yields 


-fv 

Jr 




G r K„ G, 
. —K„G r 





(13.57) 


This can be put in the more compact form 



with T = [ —G,. I] and T q = f G q 0]. 


q, + 


^ q = T r K„ Tu-T 7 q, + T T q q, 


(13.58) 


The end result is that the free-free stiffness is T 7 F rr ' T. Alternatively (13.58) may be derived by plugging 
d,. = u,. — G, u v into (13.52) and then into (13.48). 

Remark 13.6. Let S be a n s x n s nonsingular matrix. The matrix R built by the prescription 


R = 



(13.59) 


is called a rigid body motion matrix or simply RBM matrix. The columns of R represent nodal values of rigid motions, 
hence the name. The scaling provided by S may be adjusted to make R simpler. The key property is T R = 0 and thus 
KR = T r K,, T R = 0. Other properties are studied in [258]. 

§13.4.5. Applications 


Stiffness Equilibrium Tests. If one injects u,. = Gu s and q ( = 0 into (13.58) the result is f,. = 0 and f s = 0. 
That is, all node forces must vanish for arbitrary u s . This test is useful at any level (element, superelement, full 
structure) to verify that a directly generated K (that is, a K constructed independently of overall equilibrium) 
is “clean” as regards rigid body modes. 

Element Stiffness from Flexibility. Here F,.,. is constructed at the supported element level, inverted to get K„. 
and rigid motions injected through (13.58). Applications to element construction are illustrated in §13.5. 

Experimental Stiffness from Flexibility. In this case F,.,. is obtained through experimental measurements on 
a supported structure or substructure. 11 To insert this as a “user defined superelement” in a DSM code, it is 
necessary to produce a stiffness matrix. This is done again by inversion and RBM injection. 


§13.5. Flexibility Based Derivations 

The equilibrium theorems of the foregoing section are applied to the flexibility derivation of several one¬ 
dimensional elements. 

§13.5.1. Timoshenko Plane Beam-Column 

A beam-column member combines axial and bending effects. A 2-node, straight beam-column has three 
DOFs at each node: the axial displacement, the transverse displacement and a rotation. If the cross section is 
doubly symmetric, axial and bending effects are decoupled. A prismatic, plane element of this kind is shown 
in Figure 13.15(a). End nodes are 1-2. The bending component is modeled as a Timoshenko beam. The 
element is subjected to the six node forces shown, and to a uniformly distributed load q 0 . To suppress rigid 
motions node 1 is fixed as shown in Figure 13.15(b), making the beam a cantilever. Following the notation of 
§ 13.4. l—§ 13.4.2, 



r u x 11 


r u x2 1 


~ w 


fx 2 


" 0 " 


O 

O 

i —H 

U* = 

i 

_i 

, d, = u, = 

i 

1_ 

, fv = 

fy i 

_ m \ _ 

, f, = 

fy2 

_ m 2 _ 

. q = 

ql 

_ 0 _ 

, G(x) = 

0 1 x 
_0 0 0_ 


(13.60) 


li 


The classical static tests on an airplane wing are performed by applying transverse forces and torques to the wing tip 
with the airplane safely on the ground. These experimental influence coefficients can be used for model validation. 
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Figure 13.15. Flexibility derivation of Timoshenko plane beam-column stiffness: (a) element and node 
forces, (b) removal of RBMs by fixing left node, (c) FBD that gives internal forces at varying x. 


Further, G, = G(l) and G q = G(l/2). The internal forces are the axial force F(x), the transverse shear 
V(x) and the bending moment M (x ). These are directly obtained from statics by doing a free-body diagram 
at distance x from the left end as illustrated in Figure 13.15(c). With the positive convention as shown we get 

N(x) = -f x 2 , V(x) = -f y2 -q 0 (l-x), M(x) = m 2 + f y2 (i ~ x) + \q 0 (l - x) 2 . (13.61) 


Useful check: clM/dx = V. Assuming a doubly symmetric section so that N and M are decoupled, the 
element TCPE functional is 


n* 



N 2 M 2 
~EA + £7 + 



dx — f ( 7 d,. 


±f r F rr f, + f r b, 


2 r 


fjd r + n 


* 

0 ’ 


3 2 n* 

in which F. .. = - 

3f, 3f, 


r t 

EA 

0 

0 


0 

0 

t \4 + 4>) 
2AEI 

l 2 

2EI 

i — qo 

e 3 ( 4 + <h) 
12EI 

0 

i 2 

2 El 

l 

El J 


£ 2 

L 6EI J 


(13.62) 


Term fl ( * is inconsequential, since it disappears on differentiation. 

Applying the TCPE principle yields F,.,. f ( . = b, — d,. This is inverted to produce the deformational stiffness 
relation f, = K,.,. d,- + q ( ., in which 


k rr = f: 


r EA 

0 

0 "I 



i 

0 

12EI 

6 El 

, (J r — q 0 b r — 

0 

f 3 (l + d>) 

l 2 { 1 + 4>) 

q 0 l/2 

o 

6EI 

£7(4+ $) 


—qol 2 /l2 


t 2 { 1 + 4>) 

7(1 + 4)) J 




(13.63) 


To use (13.58) the following transformation matrices are required: 



"-I 

0 

0 - 



-1 

0 

0" 




0 

-1 

0 



0 

1 

0 


1 

° rS ( 

1_ 

[-?]■ 

0 

1 

-t 

0 

-1 

0 

T = 

’ 

o p 
_1 

0 

0 

i/2 

0 

1 

0 

- q = 
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0 

1 

0 



0 

0 

0 




_ 0 

0 

1 _ 



_0 

0 

0_ 




(13.64) 
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ClearAll [Le, El, GAs, <t>, qO, fx2, fy2,m2] ; GAs=12*EI/ (T*Le A 2) ; 

F=fx2; V=-fy2-q0*(Le-x); M=m2+fy2*(Le-x)+(1/2)*q0*(Le-x) A 2; 

Print["check dM/dx=V: ", Simplify[D[M,x]-V] ]; 

Ucd=F A 2/(2*EA)+M A 2/(2*EI)+V A 2/(2*GAs); 

Uc=Simplify[Integrate[Ucd,[x,0,Le}]]; Print["Uc=",Uc]; 
u2=D[Uc,fx2]; v2=D[Uc,fy2]; 02=D[Uc,m2]; 

Frr={{ D[u2,fx2], D[u2,fy2], D[u2,m2]}, 

{ D[v2,fx2], D[v2,fy2], D[v2,m2]}, 

{ D[02,fx2], D[02,fy2], D[02,m2]}}; 

br={D[Uc,fx2],D[Uc,fy2],D[Uc,m2]}/.{fx2->0,fy2->0,m2->0}; Print["br=",br]; 
Frr=Simplify[Frr]; Print["Frr=",Frr//MatrixForm]; 

Krr=Simplify[Inverse[Frr]]; Print["Krr=",Krr//MatrixForm]; 
qr=Simplify[-Krr.br]; Print["qr=",qr]; 

TT={{-1,0,0},[0,-1,0},{0,-Le,-l},{1,0,0},{0,1,0},{0,0,1}}; 

T=Transpose[TT]; Simplify[Ke=TT.Krr.T]; Print["Ke=",Ke//MatrixForm]; 

GrT={{1,0,0},{0,1,0},{0,Le,l}}; Gr=Transpose[GrT]; 

GqT={{1,0,0},{0,1,0},{0,Le/2,1}}; Gq=Transpose[GqT]; 

Print["Gr=",Gr//MatrixForm," Gq=",Gq//MatrixForm]; 
qv={0,qO*Le,0}; Print["qs=",Simplify[-GrT.qr-GqT.qv]]; 


Figure 13.16. Script to derive the stiffness matrix and consistent load vector of the prismatic, 
plane Timoshenko beam element of Figure 13.15 by flexibility methods. 

Injecting the rigid body modes from (13.58), K e = T 7 K„. T and f = T 7 q r , yields 
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f e =q 0 l[0 1/2 1/12 0 1/2 -1/ 12] r . 

(13.65) 

The bending component is the same stiffness found previously in §13.2.6; compare with (13.18). The node 
force vector is the same as the consistent one constructed in an Exercise. A useful verification technique is to 
support the beam element at end 2 and recompute K' and E. This should reproduce (13.65). 

All of the foregoing computations were carried out by the Mathematica script shown in Figure 13.16. 

§13.5.2. Plane Circular Arch in Local System 

In this and next subsection, the flexibility method is used to construct the stiffness matrix of a curved, prismatic, 
plane beam-column element with circular profile, pictured in Figure 13.17(a). The local system {.r, y } is defined 
as shown there: x is a “chord axis” that passes through end nodes 1-2, and goes from 1 to 2. Axis y is placed 
at +90° from x. No load acts between nodes. 

In a curved plane element of this nature, axial extension and bending are intrinsically coupled. Thus consid¬ 
eration of three freedoms per node is mandatory. These are the translations along x and y, and the rotation 6 
about z. This element can be applied to the analysis of plane arches and ring stiffeners (as in airplane fuselages 
and submarine pressure hulls). If the arch curvature varies along the member, it should be subdivided into 
sufficiently small elements over each of which the radius magnitude is sensibly constant. 

Care must be taken as regards sign conventions to ensure correct results when the arch convexity and node 
numbering changes. Various cases are pictured in Figure 13.18. The conventions are as follows: 

(1) The local node numbers define a positive arclength traversal along the element midline as going from 
1 to 2. The curved length i (not shown in figure) and the (chord) spanlength 2 S are always positive. 
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FIGURE 13.17. Flexibility derivation of plane circular arch element: (a) element and node forces, 
(b) removal of RBMs by fixing left node, (c) free body diagram of varying cross section. 



FIGURE 13.18. Sign conventions in derivation of circular arch element. 

(2) The arch rise H, the angular span 2 <p and the arch radius R are signed quantities as illustrated in 
Figure 13.18. The rise is the distance from chord midpoint to arch crown C: it has the sign of its 
projection on y. The angular span 2 cp is that subtended by the arch on moving from 1 to 2: it is positive 
if CCW. Finally, the radius R has the sign of cp so i = 2R(p is always positive. 

Some signed geometric relations: 


S — R sirup, H = R(costp — 1) 


(13.66) 


The location of an arch section is defined by the “tilt angle” xp measured from the circle center-to-crown 
symmetry line OC, positive CCW. See Figure 13.18. The differential arclength ds = R dxp always points in 
the positive traversal sense 1 —> 2. 

The rigid motions are removed by fixing the left end as shown in in Figure 13.17(b). The internal forces F(xp), 
V(xp) and M(xp ) at an arbitrary cross section are obtained from the FBD of Figure 13.17(c) to be 


F = fx 2 cos xp+f y2 sin xp, V = f x2 sin xp-f y2 cos xp, M = m 2 +f x2 R(cos (p— cos xp)+f y2 R(sinxp— sirup). 

(13.67) 

For typical straight beam-column members there are only two practically useful models: Bernoulli-Euler (BE) 
and Timoshenko. For curx’ed members there are many more. These range from simple corrections to BE 
through theory-of-elasticity-based models. The model selected here is one of intermediate complexity. It is 
defined by the internal complementary energy functional 


U* = 


a 


(F - M/R) 2 
TEA 


+ 


M 2 

2EI 


ds = 


/:< 


(F - M/R) 2 
TEA. 


+ 


M 2 

TEI 


R dxp. 


(13.68) 


The assumptions enbodies in this formula are: (1) the shear energy density V 2 /(2G A s ) is neglected; (2) the 
cross section area A and moment of inertia I are unchanged with respect of those of the straight member. 
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These assumptions are reasonable if \R\ > 10 r, where r — + s/l / A is the radius of gyration of the cross 
section. Further corrections are treated in Exercises. To simplify the ensuing formulas it is convenient to take 


EA 


El 

4> 2 /? 2 0 2 


AEI 

A> 2 1 2 ' 


or 



with r 2 


I 

A' 


(13.69) 


This defines T as a dimensionless geometric parameter. Note that this is not an intrinsic measure of arch 
slenderness, because it involves the element length. (In that respect it is similar to 4> of the Timoshenko 
beam element, defined by (13.17).) The necessary calculations are carried out by the Mathematica script 
of Figure 13.21, which has been pared down to essentials to save space. The deformational flexibility and 
stiffness computed are 



" F u 

Fn 

Fn' 


" K n 

K{2 

K 13 

F = 

1 rr 


F22 

F23 

, K„ = F- 1 = 


K 2 2 

K23 


_ symrn 


F33 _ 


_ symm 


K33 _ 


in which (first expression is exact value, second a Taylor series expansion in 0) 


Fn = ^|^(20(2 + 0 2 vT 2 ) + 20(l + 0 2 4/ 2 )cos20-3sin20) = (l 5 fl/ 2 + </, 2 (2 - 15fl/ 2 )) + 0(^ 4 ) 

Fi= ( sin< ft ~ ^ + <ft 2 ^ 2 )cos 4>) = r^j (l - 34> 2 ) + 0(0 3 ), 

Fl3 = iFir ( sin</> ~ 0(1 + t 2 * 2 '*™*) = §Ti (<l ~ 3 * 2 ) + 0(</,3) ’ 

Fl2 = J/03 ( 2< ^ 2 + 0 2 ^ 2 ) - 2^(1 + 0 2 fl' 2 ) cos20 - sin20) = ^ (20 + 30 2 (5fl/ 2 - 2)) + O(0 4 ), 

F 3 = 1 + (6 + 0 2 (6T /2 - 1)) + O(0 4 ), F 33 = ^j (1 + 0 2 fl/ 2 ). 


Introduce d\ = 0(1 + <p 2 A> 2 )(<p + sin0cos0) — 2sin 2 0 and d 2 = 0 — sin0cos0. Then 


(13.70) 
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(1 + 0 2 fl/ 2 ) = _ (45fl/ 2 + 0 2 ( 15fl/ 2 + 45vfl 4 - T)) + <9(0 4 ), K u = 0, 

(0(1 + 0 2 fl' 2 )cos0 - sin 0) = ^|(3fl' 2 - 1) + O(0 3 ), 

12 El 0 , fsin0 El , , 

= ^ r (5 + 0 2 )+O(0 4 ), ^ 23 -- ^ Kl2 = ~w (3O + 0 2 ) + O(0 4 ), 

(80 2 (3 + 20 2 vT 2 ) - 9 + 16cos20 - 7 cos40 - 80 sin20 (2 + (1 + 0 2 fl> 2 ) cos 2 0)) 

(l80vh 2 + 0 2 (5 - 484> 2 )) + O(0 4 ), 

(13.71) 


The constraint K 22 = —K 22 t sin 0/(20) must be verified by any arch stiffness, regardless of the TCPE form 
used. The necessary transformation matrix to inject the rigid body modes 
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(13.72) 
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Applying the congruent transformation gives the free-free stiffness 
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The only new entry not in (13.71) is 

El 

*36 = -*33 - K 23 £smcj)/cj) = — (90vh 2 + 0 2 (12i/r 2 - 5) + 0(</> 4 ) 


(13.73) 


(13.74) 


As a check, if </> 0 the entries reduce to that of the beam-column modeled with Bernoulli-Euler. For 

example *„ -» 4 EI^/l 3 = EA/l, K 22 12EI/1 3 , * 33 -> 4 El/i, K 36 2EI/Z, etc. 

§13.5.3. Plane Circular Arch in Global System 

To use the circular arch element in a 2D finite element program it is necessary to specify its geometry in the 
{.x, y} plane and then to transform the stiffness (13.73) to global coordinates. The first requirement can be 
handled by providing the coordinates of three nodes: {.x, , v,}, i = 1, 2, 3. Node 3 (see Figure) is a geometric 
node that serves to define the element mean curvature but has no associated freedoms. 

(Section to be completed). 


13-25 



Chapter 13: ADVANCED ONE-DIMENSIONAL ELEMENTS 


ClearAll [(]), l|/, T, T, R, F, M, V, Le, EA, El ] ; 

EA=4*EI/(T'2*Le A 2); 

V=-fy2*Cos [\|/] +fx2*Sin [Vj/] ; F=fx2*Cos [\|/] +fy2*Sin [V|/] ; 

M=m2+fx2*R* (Cos [\jr] -Cos [(|) ] ) +fy2*R* (Sin [\|/] +Sin [({> ] ) ; 

Ucd=(F-M/R) A 2/(2*EA)+ M A 2/(2*EI); 

Uc=Simplify [Integrate [Ucd*R, {V)/, -(|>, <|> } ] ] ; Print [ "Uc=" , Uc] ; 
u2=D[Uc,fx2]; v2=D[Uc,fy2]; 02=D[Uc,m2]; 

Frr=Simplify[{{ D[u2,fx2], D[u2,fy2], D[u2,m2]}, 

{ D[v2,fx2], D[v2,fy2], D[v2,m2]}, 

{ D[02,fx2], D[02,fy2], D[02,m2]}}]; 

Frr=FullSimplify [Frr/ . [R->Le/ (2*cf)) } ] ; Print [ "Frr=" , Frr//MatrixForm] ; 
Krr=FullSimplify[Inverse[Frr]]; Print["Krr=",Krr//MatrixForm]; 

TT= {{—1,0,0}, [0,—1,0}, [0, —2*R*Sin [<|> ] , 1}, [1,0,0], (0,1,0], [0,0,1}}; 
TT=TT/. [R->Le/ (2*4>) } ; T=Transpose [TT] ; 

Ke=Simplify[TT.Krr.T]; Print["Ke=",Ke//MatrixForm]; 


FIGURE 13.19. Script to produce circular arch element stiffness in local coordinates. 




FIGURE 13.20. Plane circular arch element in global coordinates: (a) geometric id (b) intrinsic 

geometry recovery. 


Figure 13.21. Module to produce plane circular arch element stiffness in global coordinates. 

(To be done) 
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Two-element 
patch ijk 



Trial functions 
at node j 


FIGURE 13.22. Repeating beam lattice for accuracy analysis. 


§13.6. * Accuracy Analysis 

This section presents the accuracy analysis of a repeating lattice of beam elements, analogous to that done for 
the bar element in §12.5. The analysis uses the method of modified differential equations (MoDE) mentioned 
in the Notes and Bibliography of Chapter 12. It is performed using a repeated differentiations scheme similar 
to that used in solving Exercise 12.8. Only the case of the Bernoulli-Euler model is worked out in detail. 

§13.6.1. * Accuracy of Bernoulli-Euler Beam Element 

Consider a lattice of repeating two-node, prismatic, plane Bernoulli-Euler beam elements of rigidity El and 
length i, as illustrated in Figure 13.22. The system is subject to an arbitrary lateral load q(x), which is assumed 
infinitely differentiable in x. From the lattice extract a patch of two elements: (1) and (2), connecting nodes 
i-j and j-k, respectively, as shown in Figure 13.22. 

The FEM patch equations at node j are 


£7 [-12 —61 24 0 -12 

£3 L 61 2l 2 0 U 2 -61 


Expand v (x ), 9 (x ) and q (x) in Taylor series about x = Xj , truncating at n +1 terms for v and 9 , and m +1 terms 
for q. Using xjs = (x — X;)/l, the series are v(x) — v< + xlsl v' ; + Nr 2 l 2 12\)v" ; + ... + {xj/ n i n / n\)v^, 8(x ) = 
e j + xlfie , i + (xlf 2 l 2 /2\)9" + .. . + W n t n /n\)6 ] j l \aadq(x) = q i +^lq' i + (^ 2 t 2 /2\)q'! + .. , + /m\)q [ ' n] . 

Here v l "\ etc., is an abbreviation for d" v(Xj)/dx", 12 Evaluate the v(x) and 8(x) series at i and k by setting 
xj/ = ±1, and insert in (13.75). Use the q(x) series evaluated over elements (1) and (2), to compute the 
consistent forces and m j as 



(13.76) 

Here N 3 (1) = 1(1 + § (1) ) 2 (2 - £ (1) ), = -f (1 + ? (1) ) 2 (1 - ? (1) ) n[ 2) = ±(1 - § (2, ) 2 (2 + ^ (2) ), and 

Nj 2 ’ = |(1 — ^ (2) ) 2 (1 + £ (2) ) are the Hermitian shape functions components of the j node trial function, 
whereas q (1) = q(xjr a) ), \jf a) = —1(1 — § (1> ) and g (2) = q(xj/ (2) ), i/^ (2) = 1(1 + £ (2) ) denote the lateral loads. 


12 Brackets are used instead of parentheses to avoid confusion with element superscripts. If derivatives are indexed by 
primes or roman numerals the brackets are omitted. 
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To show the resulting system in compact matrix form it is convenient to collect the derivatives at node j into 
vectors: 


V/ = [m V. < 0 , = [9, 9' 9" 


■0) n] f , q j = [qj q 


n [«] -i T 

j ‘h •••?} ] 


(13.77) 


The resulting differential system can be compactly written 
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Here S vv , S ve , Sg v and S,„, are triangular Toeplitz matrices of order (n+1) x («+1) whereas P t , and P s are 
generally rectangular matrices of order (n+ 1) x (m+1). Here is the expression of these matrices for n = 8, 
m = 4: 
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Elimination of 9j gives the condensed system 


Sv,. = Pq y , with S = S BB -SJS e _ X> P = P.-SAX' 


(13.81) 
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The nontrivial differential relations 13 given by S \j = Pqare EI(v'? — l 4 v l j"' / 720) = q, — l 4 q'? / 720, 
EIv? = q'j, EIv? 1 = q'j, EIv?" = q'”, and EIv' = q" ! . The first one is a truncation of the infinite order 
MoDE. Elimination of all vj derivatives but v l ? yields 


EIv'? = qj, (13.83) 

exactly. That this is not a fluke can be confirmed by increasing n and while taking m — n — 4. The first 4 
columns and last 4 rows of S, which are always zero, are removed to get S. The last 4 rows of P, which are 
also zero, are removed to get P. With m = n —4 both matrices are of order n —3 x n— 3. Then S = P for any 
n, which leads immediately to (13.83). This was confirmed with Mathematica running n up to 24. 

The foregoing analysis shows that the BE cubic element is nodally exact for any smooth load over a repeating 
lattice if consistent load computation is used. Exercise 13.18 verifies that the property is lost if elements are 
not identical. Numerical experiments confirm these conclusions. The Laplace transform method only works 
part way: it gives a different infinite order MoDE but recovers (13.83) as n —> oo. 

These accuracy properties are not widely known. If all beam elements are prismatic and subjected only to 
point loads at nodes, overall exactness follows from a theorem by Tong [821], which is not surprising since 
the exact solution is contained in the FEM approximation. For general distributed loads the widespread belief 
is that the cubic element incurs 0(£ 4 ) errors. The first study of this nature by Waltz et. al. [858] gave the 
modified differential equation (MoDE) for a uniform load q as 


l 4 

720 


_Q_ 
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The above terms are correct. In fact a more complete expression, obtained in this study, is 


(13.84) 


El 


v' v - v mu + 

720 


l 6 

- v 

3024 


It 


259200 




q iv + 


t 

- Q 

3024 1 


It 


259200 


+ ... 
(13.85) 


But the conclusion that “the principal error term” is of order t 4 [858,p. 1009] is incorrect. The misinterpretation 
is due to (13.85) being an ODE of infinite order. Truncation is fine if followed by elimination of higher 
derivatives. If this is done, the finite order MoDE (13.83) emerges regardless of where (13.85) is truncated; an 
obvious clue being the repetition of coefficients in both sides. The moral is that conclusions based on infinite 
order ODEs should be viewed with caution, unless corroborated by independent means. 


13 Derivatives of order 4 and higher are indicated by Roman numeral superscripts instead of primes. 
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§13.6.2. * Accuracy of Timoshenko Beam Element 

Following the same procedure it can be shown that the infinite order MoDE for a Timoshenko beam element 
repeating lattice with the stiffness (13.18) and consistent node forces is 


, ® ^ 4 (l + 50 - 5<D 2 ) f 6 (20 + 7<J>-70cb 2 + 35<I> 3 ) r 

El v ; +- v)’ --- -v v u + —- -v x , + ... 

1 ' 12 7 720 1 60480 1 


= c li 


i 4 (\ + 54>) - £ 6 (20 + 144> — 35O 2 ) . , 

70Q V + — -- -qT) + ..., in which 4> = 12 EI/(GA S £ 2 ). 


,J 60480 Jj 

(13.86) 

This reduces to (13.85) if <t> = 0. Elimination of higher order derivatives gives the finite order MoDE (aka 
FOMoDE) 

£ 2 $ „ El „ 

-q'. (13-87) 


EIv l ] v = q j -—q j =q j 


GA, 


which repeats for any n > 8. This happens to be the exact governing differential equation for a statically 
loaded Timoshenko beam [282, p. 23]. Consequently the Timoshenko beam is nodally exact under the same 
conditions previously stated for the Bernoulii-Euler model. 


Notes and Bibliography 

The material of this Chapter interweaves the very old and the very new. Before energy methods came to 
the FEM forefront by 1960 (see historical sketch in Appendix O), ordinary differential equations (ODE) and 
flexibility methods were essential part of the toolbox repertoire of the structural engineer. Traces of that 
dominance may be found in the books by Przemieniecki [665], Pestel and Leckie [637], and the survey by 
Gallagher [311]. Energy derivations were popularized by Archer [24,25], Martin [516], and Melosh [534,535]. 
For one-dimensional elements, however, results are often identical. This can provide a valuable crosscheck. 

The Legendre interpolation (13.2) was introduced in [254] to study optimal mass-stiffness combinations for 
beam elements in the context of finite element templates [246,261]. The diagonal covariance matrices Q„ 
given in (13.4) play a key role in model customization. The hinged element stiffness (13.14) is rederived in 
the Advanced FEM Lecture Notes [271]. using a mixed variational principle. The separation of uncoupled 
rigidity effects in stiffness forms such as (13.8) and (13.9) is suggested by template theory [261]. 

The Timoshenko beam model was originally proposed in [809]. Timoshenko cleverly packaged the model 
with miscellaneous ingredients introduced earlier by Bresse, Rayleigh and Hencky. It has become important 
as a tool for transient response and control simulations because its dynamic form is strictly hyperbolic. 14 
The Timoshenko beam element stiffness (13.18) first appeared in [834] in the guise of a spar element for 
use in aircraft structures; the end node freedoms of that element differring from the classical set used here. 
The particular form (13.18) is derived in Section 5.6 of [665] using ODEs. If transverse displacemenets and 
rotations are independently assumed, this beam model pertains to the class of “C° elements” that have been 
extensively studied in the FEM literature after 1968 in conjunction with the isoparametric approach. The book 
of Hughes [424] provides a comprehensive treatment of such methods for beams and plates. 

The classical work on beams on elastic foundations is by Hetenyi [393]. Useful solutions are tabulated in 
Roark-Young’s handbook [713]. 

That the use of homogeneous solutions of governing differential equations yields nodally-exact stiffness 
equations was first proven generally by Tong [821] in a Galerkin context. This derivation procedure, however, 
was rarely used after the 1960s. Two obstacles: (1) it is largely restricted to either one dimensional elements, 


14 The Bernoulii-Euler beam dynamic model is parabolic and thus exhibits an infinite transverse wave speed. Such a model 
is unsuitable for wave propagation problems. 
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or to problems with special symmetries that can be modeled with ODEs; 15 (2) rapidly increasing solution 
complexity in complicated problems discouraged hand derivations. Whereas the first limitation still holds, the 
increasing availability of CAS allows timely consideration of more difficult problems. The construction of the 
exact beam-on-Winkler-foundation element in §13.3 offers a case in point. Using Mathematica the complete 
derivation, checking and fully-symbolic testing took about 6 hours, whereas a hand derivation, coding and 
numerical testing would likely take weeks or months. The main application of exact elements appears to be 
a priori error estimation: how many simpler elements are needed to do the job of an exact one? 

The construction of stiffness matrices from flexibility information was historically one of the first techniques 
by which stiffness equations of MoM members were derived. The rigid body injection method of §13.4.4 
largely follows Section 6.6 of [665]. The presentation of discrete-system equilibrium theorems in §13.4.2 
includes a new ingredient missing from previous work: handling non-self-equilibrated loading systems. This 
extension removes the 40-year-old objection that flexibility methods (or more generally, schemes based on the 
TCPE principle) are unable to produce equivalent or consistent node forces. 

The use of equilibrium methods for multidimensional finite elements was pioneered by Fraeijs de Veubeke in 
the 1960s and early 1970s. His obsession with solution bounding got these methods seriously stuck, however, 
because of difficulties in interelement connections that maintain system-level equilibrium, as well as avoidance 
of spurious modes. More practical extensions lead to the so-called Trefftz and equilibrium hybrid methods. 
These are presently the topic of active research 16 but require advanced variational techniques beyond the scope 
of this book. Another recent advance is the discovery of the free-free flexibility as the true dual of the free-free 
stiffness [250,258,251]. This extension relies heavily on projection operators. 

That Hermitian BE element models are nodally exact if consistent loads are used is stated in [437, Sec. 8.3] as 
the “beam theorem.” Despite the name, no proof is given; only anecdotal evidence — it was likely discovered 
by numerical experimentation. A proof based on Fourier series appears in [282]. 

As the study in §13.6 illustrates, modified equation methods in boundary value problems 17 are delicate and 
should be used with care. Their intricacies baffled Strang and Fix who, upon doing Fourier analysis of a cubic 
beam element, incorrectly stated [772, p. 171] that only one of the discrete equations — that for u(x) — is 
consistent “and the others are completely inconsistent.” The alternative is the variational approach to error 
analysis. Although more robust and forgiving, predictions are often so conservative as to be of little value. 

References 

Referenced items have been moved to Appendix R. 


15 For example, symmetrically loaded circular plates or shells of revolution. 

16 Along with discontinuous Galerkin methods, a reinvention of Fraeijs de Veubeke’s weakly diffusive models. 

17 They are more forgiving in initial value problems. 
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Homework Exercises for Chapter 13 
Advanced One-Dimensional Elements 

EXERCISE 13.1 [A: 15] Evaluate the strain and stress fields associated with the Timoshenko beam displace¬ 
ment field (13.16). 

EXERCISE 13.2 [A/C:20] Find the consistent node forces for a Timoshenko beam element under an applied 
distributed moment m(£) = b 0 + Zq£ + b 2 % 2 + bjtj, where m is z moment per unit length. Hint: use the 
integral of m(§)#(§) instead of q(%) w($) in the procedure described in §13.2.8. 

EXERCISE 13.3 [A/C:20] Construct the stiffness matrix a Timoshenko plane beam element with a hinge at 
the center. Hint: set R Bs = 0, R Ba = El and R s = 12£7/(<M 2 ) in (13.13). Verify that it has rank one by 
computing its eigenvalues. 

EXERCISE 13.4 [A/C: 15=5+10] Consider a cantilever beam of constant bending rigidity El, constant shear 
rigidity GA S and length L, which is modeled with one Timoshenko element. The beam is end loaded by a 
transverse force P. The support conditions are iq = 6, =0. 

(a) Find the end displacement v 2 and end rotation d 2 in terms of P, E, G, /, A s and L. Compare with the 
analytical values P L 3 /(3EI) + PL/(GA S ) and PL 2 /(2EI), respectively. 

(b) Why does the finite element model provides the exact answer with one element? 

EXERCISE 13.5 [A:25] (Requires math ability). Discuss what happens in (13.18) if —> oo. Is the result 
useful for a shear-only “spar” element? Hint: eliminate 6\ and d 2 by a master-slave MFC. 

EXERCISE 13.6 [A: 10] For a given number of elements N e of length £ = 2 L/N e , relate x and X in Example 
13.2. 

EXERCISE 13.7 [A/C:40] (research paper level). Derive an exact Timoshenko-beam-on-Winkler-foundation 
equation method.element using the differential equation method. 

EXERCISE 13.8 [C:30] Write Mathematica code to verify the nodal exactness conclusion of §13.6.1 using 
the repeated differentiation approach. 

EXERCISE 13.9 [C:30] As above, but using the Faplace transform. Show that this only does half the job. 
EXERCISE 13.10 [A/C:35] Find the general symbolic expression of the terms in 

EXERCISE 13.11 [A/C:40] (research paper level) Analyze nodal accuracy if the length of the beam elements 
(1 ± a)i, where 0 < a < 

EXERCISE 13.12 [A/C:35] Using Mathematica, verify the results (13.86) and (13.87) in §13.6.2. 
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